

if T2C==1
    Net = importdata('Pairwise_T2C.csv');
    Indiv = importdata('Indiv_T2C.csv');
    sch_size = importdata('sch_size_T2C.csv');
else
    Net = importdata('Pairwise_T1.csv');  
    Indiv = importdata('Indiv_T1.csv');
    sch_size = importdata('sch_size_T1.csv');   
end

Net = Net.data;
G_ij = Net(:,1);
NETMISS_raw = Net(:,2);
if T2C==0
    NETMISS_raw = ones(size(G_ij,1),2);
end

Indiv = Indiv.data;
P = Indiv(:,1);
X_raw = Indiv(:,2:7);
LMH = Indiv(:,8:13);
y1y0 = Indiv(:,14:17);
y1y0_MISS = Indiv(:,18:21);


sch_size = sch_size.data;
schools = size(sch_size,1);
n = sum(sch_size);



% Create matrices
k = 0;
I_i = sparse(0,0);
I_j = sparse(0,0);

for s=1:schools
    sch_size_s = sch_size(s);
    I_i_s = kron(speye(sch_size_s),ones(sch_size_s-1,1));

    I_j_s_raw = repmat(speye(sch_size_s),1,sch_size_s);
    I_j_s = [];

    for i=1:sch_size_s
        for j=1:sch_size_s
            if i~=j
                I_j_s = [I_j_s; I_j_s_raw(j,sch_size_s*(i-1)+1:sch_size_s*i)];          
            end
        end
    end

    dimI = size(I_i,1);
    I_i = [I_i, sparse(dimI,sch_size_s); sparse(sch_size_s*(sch_size_s-1),k), I_i_s];

    I_j = [I_j, sparse(dimI,sch_size_s); sparse(sch_size_s*(sch_size_s-1),k), I_j_s];

    k = k + sch_size_s;

end

% Create flipping matrix
Flip_ij = sparse(size(I_j,1),size(I_j,1));
sch_min = 0;
sch_min1 = 0;    
for s=1:schools
    s
    sch_size_s = sch_size(s)
    Flip_ij_s = sparse(sch_size_s^2,sch_size_s^2);
    for i=1:sch_size_s
        i;
        for j=1:sch_size_s
            Flip_ij_s(sch_size_s*(i - 1) + j, sch_size_s*(j - 1) + i) = 1;
            Flip_ij_s(sch_size_s*(j - 1) + i, sch_size_s*(i - 1) + j) = 1;
        end
    end
    Flip_ij_s_new = [];
    k = 0;
    for i=1:sch_size_s
        for j=1:sch_size_s
            k = k + 1;
            if j ~= i
                Flip_ij_s_new = [Flip_ij_s_new; Flip_ij_s(k,:)];
            end
        end
    end

    Flip_ij_s_new2 = [];
    k = 0;
    for i=1:sch_size_s
        for j=1:sch_size_s
            k = k + 1;
            if j ~= i
                Flip_ij_s_new2 = [Flip_ij_s_new2, Flip_ij_s_new(:,k)];
            end
        end
    end                

    Flip_ij(sch_min1+1:sch_min1 + sch_size_s*(sch_size_s-1),sch_min1+1:sch_min1 + sch_size_s*(sch_size_s-1)) = Flip_ij_s_new2;

    sch_min = sch_min + sch_size_s;
    sch_min1 = sch_min1 + sch_size_s*(sch_size_s-1);
end

G_ij = log(G_ij);
G_ji = Flip_ij*G_ij;






% Data Matrices
X_i = sparse(I_i*[X_raw, P, repmat(P,1,size(X_raw,2)).*X_raw]);
X_j = sparse(I_j*[X_raw, P, repmat(P,1,size(X_raw,2)).*X_raw]);
dimX = size(X_i,2)

X_sum_less_ij = [];
k = 0;
k1 = 0;
for s=1:schools
    sch_size_s = sch_size(s);

    X_sum_less_ij = [X_sum_less_ij; (sch_size_s-1)^(-1)*repmat(sum(X_i(k1+1:k1+sch_size_s,:),1),sch_size_s*(sch_size_s-1),1) - X_j(k1+1:k1+sch_size_s*(sch_size_s-1),:) - X_i(k1+1:k1+sch_size_s*(sch_size_s-1),:)];

    k = k+sch_size_s;
    k1 = k1 + sch_size_s*(sch_size_s-1);
end


I_sum = sum(I_i,1);
I_sum1 = I_sum - ones(1,n);
BlockOnes = sparse(size(X_i,1),n);
k = 0;
k1 = 0;
for s=1:schools
    sch_size_s = sch_size(s);
    BlockOnes(k1+1:k1+sch_size_s*(sch_size_s-1),k+1:k+sch_size_s) = ones(sch_size_s*(sch_size_s-1),sch_size_s);
    k1 = k1+sch_size_s*(sch_size_s-1);
    k = k + sch_size_s;
end

Trans_doti = speye(size(I_i,1)) - I_i*(I_i./repmat(I_sum,size(I_i,1),1))';



mean_Xk_less_ij = (I_i*I_i'*X_j - X_j)./repmat(I_i*I_sum',1,dimX);
X_j_mean_Xk_less_ij = X_j.*mean_Xk_less_ij;
X_i_mean_Xk_less_ij = X_i.*mean_Xk_less_ij;


X_i_doti = Trans_doti*X_i;
X_j_doti = Trans_doti*X_j;
G_ij_doti = Trans_doti*G_ij;
G_ji_doti = Trans_doti*G_ji;

% Missing data indicators

NETMISS_ij = NETMISS_raw;
NETMISS_ji = Flip_ij*NETMISS_raw;
NETMISS = min(ones(size(NETMISS_raw,1),1), NETMISS_ij + NETMISS_ji);
A_MISS = ones(n,1) - min(I_sum' - I_i'*NETMISS, ones(n,1));    

y1_MISS = y1y0_MISS;
y1_educ_MISS = y1_MISS(:,1);
y1_gender_MISS = y1_MISS(:,2);



J_greater = sparse(size(I_j,1),1);
k = 0
for s=1:schools
    s
    sch_size_s = sch_size(s)
    for i=1:sch_size_s
        for j=1:sch_size_s
            if j ~= i
                k = k + 1;
            end
            if j < i
                J_greater(k) = 1;
            end
        end
    end
end    


% LMH
LMH_educ = LMH(:,1:3);
LMH_gender = LMH(:,4:6);

L_educ = LMH_educ(:,1);
M_educ = LMH_educ(:,2);
H_educ = LMH_educ(:,3);
L_gender = LMH_gender(:,1);
M_gender = LMH_gender(:,2);
H_gender = LMH_gender(:,3);

% y1 and y0
y1 = y1y0(:,1:2);
y1_educ = y1(:,1);
y1_gender = y1(:,2);
y0 = y1y0(:,3:4);
y0_educ = y0(:,1);
y0_gender = y0(:,2);


% Construct W_trans
k1 = 0;
k2 = 0;
for s=1:schools
    k1 = k1 + sch_size(s)^2;
    k2 = k2 + sch_size(s)*(sch_size(s)-1);
end
W_trans = sparse(k1,k2);
k1 = 0;
k2 = 0;
for s=1:schools
    k3 = 0;
    W_trans_s = [];
    for i=1:sch_size(s)
        for j=1:sch_size(s)
            W_trans_temp = sparse(1,sch_size(s)*(sch_size(s)-1));
            if i==j
                W_trans_s = [W_trans_s; W_trans_temp];
            else
                k3 = k3 + 1;
                W_trans_temp(k3) = 1;
                W_trans_s = [W_trans_s; W_trans_temp];
            end
        end
    end
    
    W_trans(k1+1:k1+sch_size(s)^2,k2+1:k2+sch_size(s)*(sch_size(s)-1)) = W_trans_s;
    
    k1 = k1 + sch_size(s)^2;
    k2 = k2 + sch_size(s)*(sch_size(s)-1);
end

% Construct weights for means
W_ij = exp(G_ij) + exp(G_ji);
W_ij = W_trans*W_ij;

W_ij_MAT_rep = zeros(n,n);
k = 0;   
for s=1:schools
    W_ij_MAT_rep(k+1:k+sch_size(s),k+1:k+sch_size(s)) = reshape(W_ij(k+1:k+sch_size(s)^2,1),sch_size(s),sch_size(s));
    k = k + sch_size(s);
end    
W_ij_MAT = bsxfun(@rdivide,W_ij_MAT_rep,sum(W_ij_MAT_rep,2));


% Structures and save
if T2C==1
    Cons_T2C = struct('schools',schools,'sch_size',sch_size,'n',n,'dimX',dimX,'X_raw',X_raw,'P',P,'X_i',X_i,'X_j',X_j,'X_j_doti',X_j_doti,'I_j',I_j,'I_i',I_i,'I_sum',I_sum,'I_sum1',I_sum1,'BlockOnes',BlockOnes,'mean_Xk_less_ij',mean_Xk_less_ij,'X_j_mean_Xk_less_ij',X_j_mean_Xk_less_ij,'X_i_mean_Xk_less_ij',X_i_mean_Xk_less_ij,'Trans_doti',Trans_doti,'W_trans',W_trans,'Flip_ij',Flip_ij,'J_greater',J_greater,'L_educ',L_educ,'M_educ',M_educ,'H_educ',H_educ,'L_gender',L_gender,'M_gender',M_gender,'H_gender',H_gender,'y0_educ',y0_educ,'y0_gender',y0_gender)
    NetOut_T2C = struct('G_ij',G_ij,'W_ij_MAT',W_ij_MAT,'y1_educ',y1_educ,'y1_gender',y1_gender);
    MISS_T2C = struct('NETMISS',NETMISS,'A_MISS',A_MISS,'y1_educ_MISS',y1_educ_MISS,'y1_gender_MISS',y1_gender_MISS);
    
    save('Mats_T2C.mat','Cons_T2C','NetOut_T2C','MISS_T2C')
else
    Cons_T1 = struct('schools',schools,'sch_size',sch_size,'n',n,'dimX',dimX,'X_raw',X_raw,'P',P,'X_i',X_i,'X_j',X_j,'X_j_doti',X_j_doti,'I_j',I_j,'I_i',I_i,'I_sum',I_sum,'I_sum1',I_sum1,'BlockOnes',BlockOnes,'mean_Xk_less_ij',mean_Xk_less_ij,'X_j_mean_Xk_less_ij',X_j_mean_Xk_less_ij,'X_i_mean_Xk_less_ij',X_i_mean_Xk_less_ij,'Trans_doti',Trans_doti,'W_trans',W_trans,'Flip_ij',Flip_ij,'J_greater',J_greater,'L_educ',L_educ,'M_educ',M_educ,'H_educ',H_educ,'L_gender',L_gender,'M_gender',M_gender,'H_gender',H_gender,'y0_educ',y0_educ,'y0_gender',y0_gender)
    NetOut_T1 = struct('G_ij',G_ij,'W_ij_MAT',W_ij_MAT,'y1_educ',y1_educ,'y1_gender',y1_gender);
    MISS_T1 = struct('NETMISS',NETMISS,'A_MISS',A_MISS,'y1_educ_MISS',y1_educ_MISS,'y1_gender_MISS',y1_gender_MISS);
    
    save('Mats_T1.mat','Cons_T1','NetOut_T1','MISS_T1')
end


% Clear working variables to reduce clutter
clearvars -except Restart chains sim_tol reps_chain chains_NoA Sim_reps PE_sim_reps burnin b_gap